Length scales, collective modes, and type- 1.5 regimes in three-band superconductors 
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The recent discovery of iron pnictide superconductors has resulted in a rapidly growing interest in multiband 
models with more than two bands. In this work we specifically focus on the properties of three-band Ginzburg- 
Landau models which do not have direct counterparts in more studied two-band models. First we derive normal 
modes and characteristic length scales in the conventional U{1) three-band Ginzburg-Landau model as well as 
in its time reversal symmetry broken counterpart with U{1) x Z2 symmetry. We show that in the latter case, 
the normal modes are mixed phase/density collective excitations. A possibility of the appearance of a massless 
phase-difference mode associated with fluctuations of the phase difference is also discussed. Next we show that 
gradients of densities and phase differences can be inextricably intertwined in vortex excitations in three-band 
models. This can lead to very long-range attractive intervortex interactions and appearance of type- 1.5 regimes 
even when the intercomponent Josephson coupling is large. In some cases it also results in the formation of 
a domain-like structures in the form of a ring of suppressed density around a vortex across which one of the 
phases shifts by tt. We also show that field-induced vortices can lead to a change of broken symmetry from 
U{1) toU{l) X Z2 in the system. In the type- 1.5 regime, it results in a semi-Meissner state where the system 
has a macroscopic phase separation in domains with broken U{1) and U{1) x Z2 symmetries. 



Superconductivity with two gaps associated with different 
bands was first theoretically predicted in 1959^. However it 
was not until 42 years later, with the discovery of MgB2^ 
that it started to attract a wide interest (for a recent review 
see^). Because the condensates in two band superconductors 
are not independently conserved, the system considered in^ 
share the same broken U{1) symmetry of the ground state as 
their single-component counterparts. The interband tunneling 
results in a system that attains its free energy minimum when 
the phase difference between the condensates is either zero or 
TT. Nonetheless in 1969 it was discussed that individual phases 
of the two condensate wave functions are important degrees 
of freedom, since they give rise to a new kind of collective 
excitations. These collective excitations are associated with 
the fluctuations of the relative phase of the two superconduct- 
ing components around its ground state value: the Leggett's 
mode"^, for a recent discussion see^. A report of the observa- 
tion of the Leggett's mode in MgB2 appeared very recently^. 

Another example of new physics which can arise in two- 
band system (as compared to their single-band counterparts) is 
associated with disparity of the characteristic length scales of 
density variations. That is, a single quantum vortex in a two- 
band system, should in general produce two different cores. 
As a consequence of this, there appears a regime which was 
recently termed type- 1.5 superconductivity^. In that regime 
the two characteristic length scales of density variations ^1 
and ^2 satisfy the condition ^1 < >/2A < ^2- For a subset 
of parameters in this regime there are thermodynamically sta- 
ble vortices with non-monotonic interaction. Namely, these 
vortices exhibit interaction which is long range attractive, and 
short-range repulsive. As a consequence of the long range 
attraction between vortices the system allows an additional 
"semi-Meissner" phase associated with a macroscopic phase 
separations in domains of the Meissner state and vortex states 
see Q.gJ~^^. For a detailed introduction see^^. 

In the last three years there has been a rapidly growing 
interest in multiband superconductivity with more than two 



components. The interest was sparked by the recent discovery 
of iron pnictide superconductors^^ and subsequent discussions 
that superconductivity in these systems may be described by 
a theory with more than two relevant bands ^^'^^. It was ob- 
served that the inclusion of a third band in the theory in sev- 
eral respects leads to qualitatively different physics compared 
to two-band systems^^"^^. The new physics arises from the 
fact that the presence of three or more components can lead 
to phase frustration. It results from competition of three or 
more interband Josephson coupling terms, which cannot all 
simultaneously attain the most energetically favorable phase 
locking pattern^^"^^. This frustration leads to Time Reversal 
Symmetry Breakdown (TRSB)^^~^^ (we discuss it more quan- 
titatively below). See also^"^'^^ for a different discussion of 
possible time reversal synmietry breakdown in iron pnictides. 
Here we show that phase frustration leads to a plethora of new 
phenomena in the physics of collective excitations and the 
magnetic response of the three-band Ginzburg-Landau model. 

In the TRSB phase there are no "phase-only" Leggett's 
modes. Instead there is a different kind of collective exci- 
tations: Mixed phase-density modes. These mixed normal 
modes have quite complex structure and, they can possess 
modes with large characteristic length scales even in the case 
of strong Josephson coupling. At the transition point to the 
TRSB regime, the length scale of one of the phase difference 
modes diverges, rendering one of the modes massless (as was 
also discussed recently in a London model^^, for other recent 
discussions of Leggett's modes in connection with iron pnic- 
tides see^^'^^). Note however that if the phase transition in the 
TRSB state is first order as was argued in Ref.^^, then there 
will be no massless mode. This is in contrast to two-band 
systems where increasing interband Josephson coupling al- 
ways diminishes disparities of the density variations^' ^^. In 
particular it implies that the type- 1.5 regime is possible in 
three band superconductors even in cases of quite strong in- 
terband Josephson coupling. Moreover we show that in three- 
band systems the Semi-Meissner state can represent not only 



a macroscopic phase separation in vortex and Meissner do- 
mains but also represent a macroscopic phase separations of 
domains with different broken symmetries. 



invariant phase differences by rewriting the model as: 



I. MODEL 



The minimal GL free energy functional to model a three- 
band superconductor is 



F=i(VxA)2+ ^ i|DV^,p+a,|V^,p + ift|V^, 



2' ^ ' A^ 2' ' ' ' 2 

i=l,2,3 
i=l,2,3 j>i 



(1) 



Here D = V + ieA and i/ji = IV^ile*'^ are complex fields 
representing the superconducting components. The phase dif- 
ference between two condensates is denoted (pij = (pj — (pi . 
The magnetic flux density through the system is given by 
B = (V X A) and the magnetic energy density is B^/2. Such 
a multicomponent GL free energy can in certain cases be mi- 
croscopically derived at temperatures close but not too close 
to Tc (for a review see^^). Indeed the existence of three super- 
conducting bands is not by any means a sufficient condition 
for a system to have GL expansion like that given in Eq. (1). 
However many of the questions which we consider below in 
fact do not require the system to be in the high-temperature re- 
gion where a GL expansion like Eq. (1) could in certain cases 
be formally justified. In what follows we will however use the 
minimal GL model since it provides a convenient framework 
to discuss this physics qualitatively. In the Eq. (1) the coef- 
ficients ai change signs at some characteristic temperatures 
which are generally different for all components. Below this 
temperature a^ < and the band is active. Above it, a^ > 
and the band is passive. Passive bands can nevertheless have 
nonzero superfluid density because of the interband Joseph- 
son tunneling terms r]ij\ijji\\tljj\ cos (fij. Thus it is possible 
in this model to have only passive bands, and still nonzero 
superfluid densities due to Josephson terms. In the three com- 
ponent model Eq. (1) there are additional terms allowed by 
symmetry, e.g. bi-quadratic terms in density (for a review of 
microscopic derivation of such terms from a weak-coupling 
two-band theory see^^). However the impact of these terms on 
length scales and vortex physics in three-band model is essen- 
tially the same as in the well studied two-band case^^. Since 
their role is mostly connected with a straightforward renor- 
malization of the length scales we will not repeat this analysis 
here. Instead we will focus primarily on the Josephson cou- 
plings, which can play principally different roles in two- and 
three- band cases. 

Let us first discuss the simplest London approximation 
i.e. |?/^| = const. Then one can extract gradients of the gauge- 



F = 



-[ ^ |V^.pV(^. + e ^ iV^.pA]^ 



A^2=l,2,3 IV-*I i=i,2,3 i=l,2,3 

(V((/9i-(/;2))' 



\M^2\ ,^,,_ ,. ,,2 



+ 



l^2||V^3| 
4Ei=l,2,3lV^.P 



{\/{^2-^3))' 



|V^l||V^3 



,(V((/9l-(^3))' 



4E..i,2,3iv^.r 

i=l,2,3 j>i 

4(VxA)^ 



(2) 



The first term features the phase gradients coupled to the vec- 
tor potential: this corresponds to the total current in the sys- 
tem. The rest of the terms correspond to counter-flow of 
carriers in different bands. Since there is no charge trans- 
fer in counter-flows there is no coupling to gauge fields. 
In the limit rjij = the second, third and fourth term 
describe neutral superfluid modes with phase stiffnesses 
|'0i||'0j|/[4Ei=i 2 3 IV^^P] which were studied in detail in^^. 
When Josephson terms are present they break symmetry by 
giving preferred values to the phase differences, yet the sys- 
tem can have fluctuations near these values. After this illus- 
tration of phase fluctuations, we discuss in the following the 
fluctuations within the full Ginzburg-Landau model which in- 
volves fluctuations of both phases and densities. 

Systems with more than two Josephson-coupled bands can 
e^xhibit phase frustration. For rjij < 0, a given Josephson in- 
teraction energy term is minimal for zero phase difference (we 
then refer to the coupling as "phase-locking" ), while when 
r]ij > it is minimal for a phase difference equal to n (we 
then refer to the coupling as "phase-antilocking" ). Two com- 
ponent systems are symmetric with respect to the sign change 
r]ij -^ —rjij as the phase difference changes by a factor tt, 
for the system to recover the same interaction. However, in 
systems with more than two bands there is generally no such 
symmetry. For example if a three band system has r^ > for 
all Josephson interactions, then these terms can not be simul- 
taneously minimized, as this would correspond to all possible 
phase differences being equal to tt. 



II. GROUND STATE OF A THREE BAND 
SUPERCONDUCTOR 

The ground state values of the fields \il)i\ and (^ij of system 
Eq. (1) are found by minimizing its potential energy 

i j>i 

Minimizing the potential energy Eq. (3) can not in general be 
done analytically. Yet, some properties can be derived from 



qualitative arguments. In terms of the sign of the r^'s, there are 
four principal situations: 

Case Sign of 7^12, r^is, r]23 Ground State Phases 



Frustrated 
Frustrated 



The case 2) can result in several ground states. If |7^23| "C 
1^12 1 , 1^13 1 , then the phase differences are generally (pij = 0. 
If on the other hand |r/i2|, |^i3| ^ |^23| then (p23 = tt and 
(/?i2 is either or tt. For certain parameter values it can also 
have compromise states with cpij not being integer multiples 

of TT. 

The case 4) can give a wide range of ground states, as can 
be seen in Fig. 1. As 7^12 is scaled, ground state phases change 
continuously from (— tt, tt, 0) to the limit where one band is 
depleted and the remaining phases are (— 7r/2, 7r/2). 




Ground state phases 



a) b) c) 
- -© -©- O- ■ 



d) 
-0 




Figure 1. Ground state phases of the three components as func- 
tion of 7712 (here (/?3 = fixes the gauge). The GL parameters are 
Qi = 1, /3i = 1, ?7i3 = ?723 = 3. For intermediate values of 7712 the 
ground state exhibits discrete degeneracy (symmetry is U(l) x Z2 
rather than U{1)) since the energy is invariant under the sign change 
(P2 ^ -^2, ^2, ^ -^2,- For large 7712 we get (^2 - V^s = vr imply- 
ing that iV^a I =0 and so there is a second transition from U{1) x Z2 
to U{1) and only two bands at the point d). Here, the phases were 
computed in a system with only passive bands, though systems with 
active bands exhibit the same qualitative properties except for the 
transition to U{1) and two bands only {i.e. active bands have non- 
zero density in the ground state). 

An important property of the potential energy Eq. (3) is that 
it is invariant under complex conjugation of the fields. That 



is, the potential energy does not change if the sign of all phase 
differences is changed, (^ij -^ —(fij. Thus, if any of the phase 
differences (pij is not an integer multiple of tt, then the ground 
state posses an additional discrete Z2 degeneracy. For exam- 
ple for a system with a^ = — 1, A = 1 and r]ij = 1, two pos- 
sible ground state are given by (pi2 = 27r/3, (pi3 = — 27r/3 
or ipi2 = — 27r/3, (pi3 = 27r/3. Thus in this case, the sym- 
metry is /7(1) X Z2, as opposed to U{1). As a result, like 
any other system with Z2 degeneracy, the theory allows an 
additional set of topological excitations : domain walls inter- 
polating between the two inequivalent ground states. Under 
certain conditions the system also does allow composite topo- 
logical excitations which are bound states of closed domain 
walls and vortices^ ^. 

We are interested in determining quantitatively (i) the 
ground state densities and phase differences and (ii) the char- 
acteristic length scales at which a perturbed field recovers its 
ground state values. These quantities are derived from a per- 
turbative expansion around the ground state. Consider the fol- 
lowing expansion of the fields entering the Ginzburg-Landau 
free energy functional Eq. (1), around the ground state 



ilji = {ui + ei{r)) exp {i{(pi + ^i(r))} , 
r 



(sin6>, cosO) . 



(4) 



The ground state densities and phases are denoted Ui and ipi 
respectively. Since we are interested in vortex excitations, lets 
consider an axially symmetric configuration by requiring that 
the field fluctuations ei{r), (})i{r) and a{r), depend only on the 
radial coordinate. The expansion Eq. (4) is inserted into the 
free energy Eq. (1) which is then sorted by growing orders in 
the fluctuations, namely F = F^^^ + F^^^ + F^^) + .... The 
condensation energy is given by F^^\ 



A. Ground state 

The ground state can be represented by the vector of the 
zeroth order degrees of freedom of Eq. (4), 



7^°^ = (7^1,7/2,^3,^^1,^^2,^3)^. 



(5) 



The fluctuation amplitudes are collected in the 6 entry vector 

7^^^ = (61,62,63,01,02,03)^. (6) 

The gauge field fluctuation a decouples from the other fluctu- 
ations. The term in the GL free energy which is linear in the 
fluctuations reads 



F^'^^ = ^^2uiei{ai -^ PiU^) + ^%-(^iej -^ UjCi) cos cpij 

i j>i 

+ ^ %-^i^j (0j - 0i) sin (fij . (7) 



j>i 



where (fij denote phase differences of the ground state. 
Eq. (7) is a linear (in the fluctuations) system of 6 equations. 
Since we consider fluctuations near the ground state it has 



to be zero for any arbitrary fluctuation. Indeed, by defini- 
tion, no fluctuation can decrease the energy of the ground 
state. Positive definiteness impHes that aU the prefactors of 
the fluctuations are zero. Thus expanding Eq. (7) and coflect- 
ing the prefactors of the fluctuation amplitudes gives the sys- 
tem of 6 equations which determine the ground state vector 
{ui^U2j Usj(fij(p2j ^3)^ • The system reads expHcitly 



7 



(0) 



^2 



^13 

2 
2 

n , /Q 3 , ^13 - , ^23 



= aiUi + PiUi + ^U2 COS (pi2 
= a2U2 + /32ul + -^Ui cos (^12 



^^3 cos (^13 (8a) 

Us cos (^23 (8b) 
U2 cos (^23 (8c) 



2 "^ 2 

= -r]i2UiU2 sin (^12 - r/13^1^3 sin (^13 (8d) 

= r]i2UiU2 sin (^12 - r?23^2^3 sin (^23 (8e) 

= rjisuius sin (^13 + 7?23^2^3 sin (^23 • (8f) 

Except under very specific conditions this cannot be solved 
analytically. In this paper we aim at the most general structure 
of the ground state, so no further assumptions will be made 
and the problem is solved using numerical methods (we here 
used Newton-Raphson algorithm). For numerical calculations 
of the ground state values of the fields, it is convenient to fix 
the gauge by for example imposing ipi = 0. 



B. Length scales 

Once the ground state 7^^^ is known, relevant informa- 
tion about the physics of the system can be extracted from 
the quadratic order F*^^^ of the fluctuation expansion (note 
that this is equivalent to considering linearized GL equations). 
The fluctuations are described by a system of Klein-Gordon 
equations for the 6 condensate fluctuations (3 densities plus 
3 phases), supplemented by a Proca field equation which de- 
scribes fluctuations of the gauge field. For studying the system 
it may be convenient to switch to a slightly different basis 



7 



(1) 



(ei,e2,e3,7ri,7r2,7r3)^ where. 



Uj 



(9) 



since in this basis, the (squared) mass matrix of the Klein- 
Gordon system is symmetric. The results will be straightfor- 
wardly switched back to the basis (j). The total functional at 
this order reads 



F(2) = & 



where 



E, 



Klein-Gordon 



^(7^^)')^ 



• Ep, 



y(i)X2 (1) 



(10) 



1,0', 2 e"^ 



E^ 



(11) 



Here ' denotes the differentiation with respect to the radial 
coordinate r. The (squared) mass matrix M^ of the Klein- 



Gordon system can easily be read from 






2 V Uj Ui 



j cos (^ij I , 



(12) 



simply by identifying the prefactors of the perturbations and 
filling the corresponding entries in the mass matrix. Before 
discussing in detail this mass matrix, let us consider the Proca 
equation, for the mass of the gauge field. It is the easiest 
length scale to derive, since the Proca equation for the gauge 
field fluctuation Eq. (11) decouples from all other. The Lon- 
don penetration depth of the magnetic field A is the inverse 
mass of the gauge field, namely 



A; 



1 



1 



^ve; 



(13) 



Length scales associated with condensate's degrees of free- 
dom are obtained in a more complicated way. Indeed they 
are given by the eigenvalue spectrum of a system of 6 cou- 
pled (static) Klein-Gordon equations, whose (squared) mass 
matrix A^^ is derived from Eq. (12). It may be instructive 
to have this mass matrix explicitly. First of all, let us remark 
that fluctuations can be separated in two groups, the 'density 
amplitude' / = (/i, /2, /3)^, and the 'normalized phase am- 
plitudes' TT = (tti, 7r2, 7r3)^. This mass matrix is a real sym- 
metric matrix, which is not diagonal and whose eigenvalues 
are the (squared) masses of the normal modes. The eigen- 
spectrum of A^^, defines the (squared) masses of the physical 
modes. The inverses of each of the masses gives the charac- 
teristic length scales of the theory. For example in a single 
component theory the inverse mass of the fluctuations of the 
modulus of the order parameter | V^| is the coherence length (up 
to a factor of >/2). In a two-component theory the fluctuations 
in the phase difference (the Leggett's mode) are characterized 
by a mass, the inverse of which sets the length scale at which a 
perturbed phase difference recovers its ground state values. In 
two-component models the density modes are mixed: i.e. the 
characteristic length scales of the density fields are associated 
with the linear combinations of the fields^'^^'^^. Physically this 
means that disturbing one density field necessarily perturbs 
the other. It also implies that, say in a vortex, the long-range 
asymptotic behavior of both density fields is governed by the 
same exponent, corresponding to a mixed mode with the low- 
est mass. 

We see that in the three-component case a new situation can 
arise where different collective modes are possible which are 
associated with mixed density and phase modes In the basis 
(/, tt), the (squared) mass matrix can be written in terms of 4 
sub-matrices 







(14) 



Where Mff and M^^ are the self-coupHng of density and 



phase fluctuations, while Mfj^ and M^^f blocks control the 
mixing of density modes and phase modes. 



Mff = 






?>p2y? 



2 m 

as + SPsul 



Mf^ = M^f = 



-fjis 



r]i2 



ms 



rii2Ui—ri23U3 

U2 



^23 

^13^1+^23^2 
Us 



2 ^2^712+^3^713 
Ui 

f]l2 

-fjis 



-r]i2 

2 '^1^12+^3^23 
U2 

-^23 



-^13 
-^23 

2 ^1^13+^2^723 



(15) 



where for having more compact expression we introduce new 
notations fjij = ^ cos (pij , and fjij = ^ sin (pij . Finally in 
order to derive the length scales associated with the conden- 
sate fluctuations, one has to diagonalize the matrix M^. Its 
eigenspectrum is the set of 6 squared masses A^f , whose cor- 
responding lengths £i = I /Mi are the physical length scales 
of a three band superconductor. In appendix A we also show 
how these length scales are expressed in different units. There 
is a spontaneously broken U{1) symmetry associated with the 
simultaneous equal changes of all phases. The mass of this 
mode is zero and the eigenvector associated with this U{1) 
zero mode can easily be decoupled. Thus one can reduce the 
size of the system. However we prefer not to decouple this 
mode from the mass spectrum, since it provides a measure of 
the error of the numerical resolution of masses of other modes. 
The corresponding degree of freedom is described by the first 
term in Eq. (2), it is a U{1) Goldstone boson which, due to 
its coupling to the gauge field A yields a massive vector field 
with the mass mproca via the Anderson-Higgs mechanism. 

Unfortunately the eigenbasis of M^ cannot be known an- 
alytically, in the general case. We calculate it numerically 
below. 



C. Numerical results 

Fig. 2 shows the ground state, eigenspectrum and eigen- 
vectors of the (squared) mass matrix in a frustrated three- 
band superconductor as a function of the Josephson couplings. 
The coupling 7/12 is fixed, while the horizontal axis gives the 
coupling coefficients 7^13 and r]23- Each eigenvector is a lin- 
ear combination of the degrees of freedom that comprises a 
physical mode, which variation length scale is given by the 
square root of the inverse of the corresponding eigenvalue in 
the eigenspectrum. The system crosses over from U{1) to 
U{1) X Z2 TRSB state at 7^13 = 7^23 ^ -3.69. Interestingly, 
in the [/(I) regime, the density modes are mixed. However, as 
can be seen from the eigenvectors, there is no mixing between 
density modes and the phase modes. Thus, perturbations of 
the densities and of the phases recover independently from 
each other. The fluctuations of the phase modes are the three- 
component generalization of the standard Leggett's modes. In 



I 

the [/(I) X Z2 regime the situation is opposite, and all eigen- 
vectors are mixed in density and phase. This indicates that 
any perturbation of the densities creates a perturbation to the 
phases, and vice versa. 

There is a point where a Leggett's mode becomes massless, 
as was also pointed out recently in the phase-only model in^^. 
This occurs at the transition from /7(1) to U{1) x Z2 (note 
however, that the transition between these states can be first 
order as discussed in^^). In Fig. 2 panel (c) the eigenvalue 
5 does indeed go to zero, indicating that the mass vanishes. 
The corresponding eigenvector can be seen in panel (h). In 
the U{1) regime it corresponds to perturbation of the phases 
1 and 2. The physical implication is that the recovery of a 
perturbation at this point is governed not by an exponential, 
but by a power law. It is only a point in the parameter space 
where this mass is zero. However there is a finite area in the 
parameter space around that point where although the mode is 
massive, the length scale associated with this it is anomalously 
large as a consequence of the frustration between Josephson 
couplings. 



III. VORTEX MATTER IN THREE-BAND TYPE-1.5 
SUPERCONDUCTORS 

A. Topological defects in three-band Ginzburg-Landau model 

Lets us start with outlining the basic properties of the vor- 
tex excitations. In case of a [t/(l)]^ Ginzburg-Landau model 
(i.e. when rjij = 0) there are three "elementary" vortex exci- 
tations associated with 27r winding in only one of the phases 
: f^ V(/?i = 27r, where a is a closed path around a vortex 
core. Such a vortex carries a fraction of flux quantum as can 
be seen from the following argument^^'^^ : the supercurrent in 
case when there is a phase winding in only one phase is 

J. = f [v^rvv^. - ^^v^^:] -e^Y. i^^i'^ ^16) 

k 

Expressing A via gradients and choosing the contour a far 
from the vortex core gives the following equation for the mag- 
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Figure 2. Ground state, eigenspectrum and eigenvectors of the mass matrix. The x-axis gives the two parameters 7713 = 7723 while the other 
parameters areai = —3, /3i = 3, a2 = —3, /32 = 3, as = 2, /Ss = 0.5, 7712 = 2.25. The eigenvectors are sorted according to 
corresponding eigenvalue, starting with the largest. The smallest eigenvalue is the zero mode associated with the spontaneously broken U(l) 
symmetry. At 7713 = 7723 ~ —3.69 there is a transition from U{1) to U{1) x Z2. The eigenvalue 5 (c) becomes zero at the transition point, 
so there appears a divergent length scale at this point which correspond to the eigenvector shown on panel (h) i.e. the phase difference mode 
becomes a scaleless collective excitation. Observe that in the U{1) region, the eigenvectors exhibit no mixing between densities and phases, 
whilst in the U{1) x Z2 region there is in fact not a single eigenvector that is not a mixing of phases and densities. Then, perturbation of 
densities are generically associated with perturbations of the phase differences in this regime. 
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where <l>o is a flux quantum. Such a fractional vortex in the 
[/7(1)]^ case has logarithmically divergent energy. Thus in 
external field a bulk three-component superconductor should 



form "composite" integer flux vortices which have phase 
winding in all components: §^ V^i = 27r, ^^ ^^2 = 
27r, f^ S/ifs = 27r. When Josephson coupling is non-zero, 
then energy of a fractional vortex diverges linearly^^ and thus 
a single integer flux vortex in a bulk superconductor can be 
viewed as a strongly bound state of three co-centered frac- 
tional flux vortices. Note that such a bound state will in gen- 
eral have three different sizes of vortex cores. The characteris- 
tic length scales of the density recovery in the vortex cores are 



determined by the inverse masses of normal modes calculated 
above. Note also that the role of Josephson interaction on vor- 
tices is different in the presence of domain walls in three-band 
U{1) X Z2 superconductors. Immediately at the domain wall 
the Josephson terms have energetically unfavourable values 
of the phase differences. Thus, if a composite vortex is placed 
on such a domain wall, the Josephson interaction can force a 
splitting of this vortex into fractional flux vortices, because the 
splitting will allow to attain a more favorable configuration of 
the phase differences^^. 



B. Qualitative argument on the vortex states in frustrated 
three-band superconductors 

The ground state of a phase frustrated superconductor is in 
many cases non-trivial, with phase differences being compro- 
mises between the various interaction terms. Inserting vor- 
tices in such a system can shift the balance between different 
competing couplings, since vortices can in general have differ- 
ent effects on the different bands. In particular, since the core 
sizes of vortices are not generally the same in all bands, vortex 
matter will typically deplete some components more than oth- 
ers and thus can alter the preferred values of the phase differ- 
ence. So the minimal potential energy inside a vortex lattice 
or cluster may correspond to a different set of phase differ- 
ences than in the vortexless ground state. In this subsection 
we give a qualitative description of it, using an ansatz-based 
argument. In the following section we study this question nu- 
merically without involving an ansatz. 

The qualitative argument is as follows. Consider the phase- 
dependent potential terms in the free energy Eq. (1) which are 
of the form 



r]ijUiUjfi{T)fj{T) COs{(pij{T)) , 



(18) 



where Ui are ground state amplitudes and each fi (r) represent 
an ansatz which models how superfluid densities are modu- 
lated due to vortices. Consider now a system where N vortices 
are uniformly distributed in a domain Q. The phase dependent 
part of the free energy is 



U^ 



i>j 



UjUj 



dr/,(r)/,(r)cos((^,,(r)). (19) 



If ifij is varying slowly in comparison with the inter vortex 
distance, then it can be considered constant in a uniform dis- 
tribution of vortices (as a first approximation). In that case 
Eq. (19) can be approximated by 

(20) 
If on the other hand (^ij varies rapidly, then it is not possible to 
define fjij without a spatial dependence. Then (^ij will depend 
on fjij (r) which is related to the local modulation functions 
fifj and vary with a length scale given by the mass matrix 
Eq. (12). 



Thus, fj is the effective inter-band interaction coupling re- 
sulting from density modulation. Since in general, fi ^ fj 
(unless the two bands z, j are identical), one must take into ac- 
count the modulation functions fi when calculating the phase 
differences. In particular, if the core size in component i 
is larger than in component j, then J drfifk < J drfjfk 
and therefore the phase differences (pij minimizing Eq. (20) 
depend on fi, and consequently on the density of vortices. 
Roughly speaking, introducing vortices in the system is equiv- 
alent to relative effective decrease of some of the Josephson 
coupling constants. 

Because the problem is nonlinear, the modulation functions 
fi generally depend on (pij since the vortex core shape de- 
pends on the inter band interactions. As a result, an exact solu- 
tion to this problem can only be found by numerical methods. 
Below, we address this problem by finding numerically vortex 
clusters solutions. Some qualitative statements can nonethe- 
less be made about these systems: 

• If band i is associated with larger vortex cores than band 
j, then with increasing density of vortices, the effective 
Josephson coupling fjik is depleted faster than fjjk- 

• The average intercomponent phase difference in a vor- 
tex cluster depends on the parameters fjij . So the inter- 
component phase differences inside a vortex cluster can 
be different from the vortexless ground state. Supercon- 
ductors with U{1) X Z2 symmetry and disparity of core 
sizes will therefore generally exhibit perturbation of the 
phase differences due to vortices. 

• The symmetry of the system depends on the inter-band 
interactions, so vortex matter can induce a phase transi- 
tion between U{1) and /7(1) x Z2 states or vice versa 

This physics depends on the spatial distribution of vortices in 
the system. 

If vortices are uniformly distributed in the sample, as is 
generally the case in clean type-2 superconductors, then the 
effective inter-band interaction strengths fjij are depleted in 
the same way everywhere in the sample. A change in symme- 
try [/(I) -^ U{1) X Z2 would then occur in the whole system 
at a certain value of applied external field. 

It also opens a possibility of a type- 1.5 regime qualitatively 
different in three-band systems than in their two band counter 
parts. Indeed, because of the non-monotonic interactions, the 
superconductor possesses macroscopic separation of Meiss- 
ner domains and vortex clusters. In the three-band case, these 
phases can exhibit different broken symmetries. For example, 
Meissner domains with the /7(1) symmetry and vortex clus- 
ters having a different symmetry, U{1) x Z2. The U{1) x Z2 
broken symmetry arising here because of the renormalization 
by vortices of the effective coupling constants fjij . If there 
is a symmetry change U{1) -^ /7(1) x Z2 associated with 
vortex clusters in the system then there will be two kinds of 
vortex clusters corresponding to with Z2 states. They will co- 
exist with the Meissner states voids which do not have the 
broken Z2 symmetry. Clearly, because of this additional dis- 
crete symmetry, inter-cluster interaction should generally be 



affected by whether the clusters are in the same, or in differ- 
ent Z2 states. When the magnetic field increases vortex clus- 
ters will merge and the entire system will be in the state with 
broken U{1) ^ U{1) x Z2 symmetry. 



C. Numerical results 

We used numerical computations to examine the questions 
which were raised about vortex matter in the previous sec- 
tions. The free energy functional Eq. (1) is minimized in 
presence of vortex matter. In these simulations the variational 
problem was defined using a finite element formulation pro- 
vided by the Freefem++^^ library framework, using a Non- 
linear Conjugate Gradient method. Reader interested in more 
technical details can refer to the appendix B. From this numer- 
ical data, several observations about vortex matter in three- 
band systems can be made. 



1. Vortex clusters with broken Z2 symmetry 

We have simulated vortex clusters in a type- 1.5 regime in 
the system given in Fig. 2 for 7^13 = 7723 = —3.7, i.e. in the 
U{1) region but close to the transition to broken time rever- 
sal symmetry. Thus, if the vortex core size in component 3 
is larger than in the bands 1 and 2, then we should expect the 
breakdown of time reversal symmetry, for a sufficiently high 
density of vortices. Fig. 3 shows that this is indeed the case. 
In the ground state, all phases are equal ((pi = (p2 = ^3), but 
once vortex clusters are present, these phases are no longer 
preferable and two other equivalent phase locking states de- 
velops. As the density of the third band is depleted, phase 
differences come to be increasingly dominated by the inter- 
band coupling between the two other bands. This coupling 
term not being minimal for 991 = (p2. 



2. Long-range intervortex forces 

Vortex matter in this system is associated with substantial 
variations of the intercomponent phase differences. As dis- 
cussed above, in three-band system there is a phase difference 
mode that becomes less and less massive as we approach the 
transition to a TRSB state. Thus in the vicinity of this point 
the mass of the corresponding mode can be very small and 
then characteristic lengths of its variation, very large. This 
provides an additional mechanism that can lead to vortex in- 
teractions at very large distances. Fig. 4 displays the same sys- 
tem as in Fig. 3, but with two vortex clusters rather than one. 
A clearly visible perturbation of the phase differences extends 
from the clusters well outside the region with magnetic field 
and far beyond the area with significant density suppression, 
providing a mechanism for long range inter cluster interaction. 



3. Vortex fractionalization in clusters 

Fig. 3 and Fig. 4 also exhibit flux fractionalization. As pre- 
viously mentioned, the model Eq. (1) allows fractional vortex 
solutions, where only one of the phases (^i winds 27r around 
some point while the rest doesn't. The flux carried by a single 
fractional vortex is given by Eq. (17). Two forces hold frac- 
tional vortices together as a one-quantum composite vortex, 
in the three-component model. First is the interaction with 
the gauge field, which gives logarithmic interaction at long 
range^^'^^. The second is the Josephson coupling, which is 
asymptotically linear. In non-frustrated superconductors the 
Josephson coupling gives attractive interaction between frac- 
tional vortices, but in frustrated systems this interaction can 
be repulsive, resulting in fractionalization of vortices^ ^ 

Consider the system in Fig. 3 and Fig. 4. The ground state 
corresponds to (pi = (^2 = ^3. Since there is an energy cost 
associated with gradients of the phase difference, these are 
expected to change slowly. Thus, far away from the cluster, 
the state is simply the ground state. Deep inside the cluster 
phase differences attain a broken U{1) x Z2 state, depend- 
ing on the density of vortex matter. If the vortex density is 
very high, then |?/^3| is very small, and we expect inside the 
cluster (^12 ^ 'T^ (provided that the cluster is large). While 
(/P12 varies slowly, the density in \^l)^\ recovers more rapidly 
at the boundary of the cluster. Thus, there may be an area 
where |(/?i2| < 7r/2 while |?/^3| is small. Consequently the in- 
teraction between fractional vortices in the bands 1 and 2 due 
to Josephson coupling is repulsive in this area. Also when 
the magnitude \iIj^\ is very small or zero, the Josephson inter- 
band coupling ?/^23 and 7/^13 which provides attractive interac- 
tion between the fractional vortices is weaker or essentially 
disappears. Thus, the interaction of the fractional vortices 
is governed by the coupling to the gauge field, which gives 
attractive interaction, and the remaining Josephson coupling, 
which in this case gives repulsive interaction. As a result, in 
that region the integer flux vortices split into fractional ones. 

This effect is found in numerical simulations of vortex clus- 
ters. Looking carefully at Fig. 3 and Fig. 4 we can see that 
the vortex cores in the bands 1 and 2 do not generally coin- 
cide. From panel g) we can read off that, at the boundary 
of the cluster, the phase difference between the components 
1 and 2 is given by |V^i||V^2| sin((/:?i2) ~ -0.7 -^ \ipi2\ ^ 
0.5 -^ cos((/:?i2) ~ a/3/2 > 0. Thus, in that region, the 
Josephson term associated with components 1 and 2 gives a 
positive energy contribution resulting in repulsive interaction 
between fractional vortices in components 1, 2 leading to frac- 
tionalization of vortices. Indeed, fractionalization occurs for 
all vortices expect those in the center of the large 8 or 9 quanta 
clusters. We observe in large systems, that fractionalization is 
important at the boundary of the clusters and becomes less 
pronounced for vortices located deep inside. The magnetic 
field is significantly smeared out, as a result of this fractional- 
ization. 

The fractionalization at the cluster's boundary has a similar 
origin as the physics which stabilizes topological solitons in 
the TRSB states in three-band superconductors^^. The differ- 
ence is however that the topological solitons discussed in^^ are 
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Figure 3. Vortex cluster exhibiting an internal Z2 state in a frustrated three-band superconductor. The panel displays the following quantities, 
a) is the magnetic field, and the densities of the different condensates are b) |^i |^, c) |^2 1^, d) |^3 1^. To monitor the relative phase differences, 
we use e) {ipi \ \ip2 \ sin((^i2), f) iV^i | iV^s | sin((/?i3), g) 1-02 1 \ip3 \ sm{(p23)- Panel h) is the energy density and the supercurrents in each condensate 
are displayed on i) Ji, j) J2 , k) J3. The GL parameters used for this simulation areai = — 3, /3i=3, a2 = — 3, /32 = 3, 0^3 = 2, /33 = 
0.5, ?7i2 = 2.25, ?7i3 = —3.7, 7723 = —3.7. Thus, they correspond to the U{1) region in Fig. 2, but close to the transition point to U{1) x Z2 
symmetry. In the ground state, all the phases are locked ((fi — (p2 — ^3) as a consequence of the Josephson couplings 7712 = 7713 = —3.7 
dominating the interaction. Inside the vortex cluster the third condensate is depleted, so the coupling terms ?7i3|^i| |V^3| cos((/?i3), {i = 1, 2} 
become much weaker while the term | V^i 1 1^2 17712 cos((/?i2) becomes dominant. In sufficiently dense vortex matter, the ground state is changed 
due to the dominating antilocking interaction between the components 1 and 2. This results in a [/(I) x Z2 state inside the vortex cluster, 
as can be seen from the phase differences plots shown in panels e), f),g). ( Note that in the very center of the vortex cluster this quantity is 
small because of small values of the prefactors |^i | |V^j |.) A closer inspection of the density panels b) and c) reveals that vortex cores in both 
densities do not necessarily superimpose (it can also be seen from the supercurrents on panels i and j) and so they are fractional vortices. This 
fractionalization occurs at the boundary of the cluster, while the vortex in the middle is a composite one-quantum vortex. 



stable bound states of Z2 domain walls and fractional vortices, 
while here there is not a Z2 domain wall, but fractionalization 
comes as a result of complicated behaviour of the fields at 



a cluster boundary which is an interface between U{1) and 
U{1) X Z2 states. 
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Figure 4. Interacting vortex clusters with internal Z2 symmetry in a frustrated superconductor. The snapshot represents a state where the 
energy is well minimized with respect to all variables except the relative positions of the weakly interacting well-separated clusters. The GL 
parameters and displayed quantities in the panels are the same as in Fig. 3. The analysis of the eigenvalues in Fig. 2 shows that there a mode 
with a very small mass, associated with the eigenvector [0, 0, 0, 1, —1, 0]. It corresponds to the mode which associated with phase difference 
fluctuations and it has the largest recovery length scale. This is indeed visible in the plots e), f), and g). The phase difference (pi2 (e) recovers 
much more slowly than the magnetic field (a) and the condensate densities (b-d). Clusters clearly interact at a distance greatly exceeding the 
length scales of density modulation and the magnetic penetration depth, as this mode stretches out between them. 



4. n -walls 



Another phenomena associated with frustrated supercon- 
ductors are objects which we term "7r-walls". In certain pa- 



by a domain- wall-like object with substantially suppressed 
density across which the phase of one of the condensates 
jumps by tt. 

An example of such an object is displayed in Fig. 5. The 



rameter regions, vortices and vortex clusters are surrounded density in the third band is small in comparison to the other 
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Figure 5. A vortex cluster surrounded by a 7r-wall. Here again displayed quantities and panel labels are the same as in Fig. 3. The GL 
parameters are ai = — 1, /3i = 1, a2 = 1, A = 0.5, as = 3, /^s = 0.5, 7712 = —2, 7713 = 2.7, 7723 = —4. In the ground state, the 
phases are locked ((fi — ^^i — (^3), but frustration occurs as 7713 = 2.7 gives an antilocking interaction (i.e. the term \i\)\ \ \\\)s \r\\s cos((/?i3) is 
minimal for (^13 = tt). As vortices are introduced in the system, the superfluid densities are depleted. It is clear from visual inspection that the 
vortices in the second band are larger than those of the first. Thus, the effective coupling 7)23 decreases faster than f]\z and so inside the vortex 
cluster the preferred phase becomes ^\ — ^2 — ^?>^ tt. Since the third band has much smaller density than the other bands, the energetically 
cheapest way of coping with this is to create a domain wall-like object where {ipsl becomes very small. It does not cost much energy to have 
a large phase gradient density, so that ip3 quickly picks up a 7r-shift in its phase. As a result the density of |^3 1 is suppressed not only in the 
vortex cores but also in a ring surrounding the vortex cluster as can be clearly seen in panel (d). 



bands. The Josephson coupling 7^12 = —2 results in locked 
phases (pi2 = 0. The system is frustrated, since 7^23 = — 4, 
preferring phase locking with respect to (/?23, and 7/13 = 2.7 
preferring phase antilocking with respect to (/913. When there 
are no vortices in the system, the term |V^2 llV^s 1^23 cos (/:?23 
dominates over |V^i||V^3|/?i3 cos(/?i3, and the ground state is 
^1 = ^2 = ^3. However, when vortices are present in the 



system, this is not necessarily the case. The vortex cores in 
the second band are larger than those of the first, and conse- 
quently, the effective coupling strength 7)23 is diminished at a 
higher pace than 7)13. Thus, inside a vortex cluster, the po- 
tential energy is minimal when (pis = (p23 = tt. To comply 
with these requirements, the system forms a domain wall-like 
object where \ijjs\ goes close to zero, and (p picks up an extra 
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phase of tt. For this particular set of parameters, this in fact 
happens even for a single vortex, as can be seen in Fig. 6. 



IV. CONCLUSIONS 

Recently there has been a growing interest to three-band su- 
perconductivity sparked by the discovery of the iron pnictide 
superconductors. The precise information about the character- 
istic parameters for these materials is not know yet. Also the 
current experiments suggest that physics of vortex ordering 
patterns in currently available samples is substantially affected 
by strong pinning^"^"^^. We presented here a general study 
showing that in a three-band system there are many phenom- 
ena which are not present in two-band models. As was previ- 
ously observed^^"^^'^^, in the presence of more than two bands, 
a system can exhibit frustration between different competing 
interband Josephson terms. We considered possible physical 
consequences using a three-band Ginzburg-Landau model. To 
observe this physics in experiment in fact does not necessary 
require a three-band superconductor but it would be sufficient 
to have a superconductor with phase anti-locking Josephson 
interaction (i.e. r] > 0). Then as was observed in^^'^^ a phase- 
frustrated state can be induced in a Josephson coupled bi-layer 
made of this and singe-band superconductors. In that case of 
the Josephson couplings is just real-space interlayer coupling. 
Thus it provides an opportunity to tune its value. 

We discussed that this can result in the appearance of modes 
with very long characteristic length scales even when the in- 
terband Josephson coupling is strong. Here we also discussed 
that in the TRSB [/(I) x Z2 state of the three-band Ginzburg- 
Landau model there are no "phase-only" Leggett's modes, but 
instead the system has different mixed phase-density collec- 
tive modes which involve both phase and density fluctuations. 
The physics of the coupled modes and associated different 
length scales substantially affects vortex matter in the system. 
The vortices can interact at distances much larger than the 
length scale of magnetic field localization or length scale at 
which most of the condensate density is recovered, because of 
the existence of slowly varying phase difference and low-mass 
mixed density modes. This can give rise to non-monotonic in- 
tervortex interaction and type- 1.5 regimes in systems where 
it would not be expected. In particular, if a IsLVge-hz parame- 
ter is estimated from the second critical field of the system, 
then this does not prohibit the existence of modes with length 
scales that substantially exceed the penetration depth even at 
strong Josephson coupling. 

Moreover the competing interactions can qualitatively af- 
fect the vortex structure as well. We showed the existence of 
vortex solution where density is suppressed not only in the 
core but also takes a second dip in some belt-like area around 
the vortex core or around the vortex cluster. Such features can 
in principle be detected in STM measurements. 

Furthermore we showed that subjecting a three band sys- 
tem to an external field which induces vortices can shift the 
balance in competing interactions and result in change of the 
ground state symmetry. In type-2 systems where vortices are 
uniformly distributed, changes in the phase difference will 



also occur quite uniformly there could be a phase transitions 
between U{1) and U{1) x Z2 states resulting from an applied 
magnetic field. In the case of type- 1.5 superconductivity, sys- 
tems will feature not only macroscopic phase separation be- 
tween vortex clusters and domains of Meissner state, but also 
a macroscopic phase separation between the domains ofU{l) 
and U{1) X Z2 ground states. The transition from the semi- 
Meissner to vortex states in that case will then be associated 
with change of the symmetry from U{1) to U{1) x Z2. 

This work is supported by the Swedish Research Council, 
and by the Knut and Alice Wallenberg Foundation through the 
Royal Swedish Academy of Sciences fellowship and by NSF 
CAREER Award No. DMR-0955902. 
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Figure 6. A single vortex in a system exhibiting vr-wall solutions. The interband coupling coefficients are 7712 = —2, 7713 = 2.7, 7723 = —4. 
Displayed quantities are here the magnetic flux a) and densities of each condensate b) |^i |^, c) |^2 1^, d) |^3 1^. Panel e) shows the total energy 
density while f), g) and h) are the supercurrents ( Ji, J2 and J3 respectively). The panel i) shows the phase difference Lpi^. The parameters of 
the system are identical to the one shown on Fig. 5. The vr-wall can be seen from the double dip in density of the third band as can be seen in 
j), as well as from the phase difference plotted in i). Thus, ^3 is zero in the center, it recovers slightly and then drops close to zero again on a 
circular area at certain distance from the vortex center. At the second drop the phase Lp2> picks up an extra phase tt as can be seen from the plot 
of (/?i3on the panel i). 
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Appendix A: Unit system 

The Ginzburg-Landau free energy Eq. (1), is written after 
suitable rescaling. Below we give details of these rescalings, 
in order to define the various quantities in the usual dimen- 
sionful theory. In the following, let us denotes the usual di- 
mensionful quantities with accentuated fonts. Consider the 
following 
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where c is the speed of light and h the reduced Planck con- 
stant, then converting the free energy Eq. (1) to 
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Here D = V — z^ A and m is the mass of the cooper pairs 

This rescaling is also applied to the perturbative expansion 
Eq. (4) of the problem, so that the Klein-Gordon system be- 
comes 

^K,e.„-o„.„„^ 1^(7(^)0' +7'')X'7'') (A3) 

Zm 

and then, (dimensionful) length scales of the massive modes 
of the condensate are 



(A4) 



In the U{1) X Z2 regime, since all the mode are mixed, 
the length scales ^i then are related to inverse masses of the 
modes. 

London penetration depth is defined through the Proca 
equation of the gauge field, which reads now in the dimen- 
sionful system 




En, 






(A5) 



London penetration depth, which gives the exponential decre- 
ment of the magnetic field in the superconductor the reads 



component Ginzburg-Landau models. They are constructed 
by minimizing the free energy Eq. (1), from an appropri- 
ate initial guess carrying several flux quanta. We consider 
the two-dimensional problem J^ = J^ Fdx'^ defined on the 
bounded domain 1] C R^, supplemented by 'open' boundary 
conditions on dQ. This 'open constraint' is a particular Neu- 
mann boundary condition, such that the normal derivative of 
the fields on the boundary are zero. These boundary condi- 
tions in fact are a very weak constraint. For this problem one 
could also apply Robin boundary conditions on dft, so that 
the fields satisfy linear asymptotic behavior (exponential lo- 
calization) Eq. (11). However, we choose to apply the 'open' 
boundary conditions which are less constraining for the prob- 
lem in question. 'Open' boundary conditions also imply that 
vortices can easily escape from the numerical grid, since it 
would further minimize the energy. To prevent this, the nu- 
merical grid is chosen to be large enough so that the attractive 
interaction with the boundaries is negligible. The size of the 
domain is then much larger than the typical interaction length 
scales. Thus in this method one has to use large numerical 
grids, which is computationally demanding. At the same time 
the advantage is that it is guaranteed that obtained solutions 
are not boundary pressure artifacts. 

The variational problem is defined for numerical compu- 
tation using a finite element formulation provided by the 
Freefem++ library^^ . Discretization within finite element for- 
mulation is done via a (homogeneous) triangulation over ft, 
based on Delaunay-Voronoi algorithm. Functions are decom- 
posed on a continuous piecewise quadratic basis on each tri- 
angle. The accuracy of such method is controlled through the 
number of triangles, (we typically used 3 ~ 6 x 10^), the order 
of expansion of the basis on each triangle (P2 elements being 
2nd order polynomial basis on each triangle), and also the or- 
der of the quadrature formula for the integral on the triangles. 

Once the problem is mathematically well defined, a numeri- 
cal optimization algorithm is used to solve the variational non- 
linear problem (i.e. to find the minima of J^). We used here 
a Nonlinear Conjugate Gradient method. The algorithm is it- 
erated until relative variation of the norm of the gradient of 
the functional T with respect to all degrees of freedom is less 
than 10-^ 



Initial guess 

The initial field configuration carrying TV flux quanta is 
prepared by using an ansatz which imposes phase windings 
around spatially separated N vortex cores in each condensates 
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Appendix B: Numerical Methods - Finite element energy 
minimization 



We provide here a detailed description of the numerical 
methods which are used to construct vortex solutions in three 



N. 



;G+iAi2 



V^3 = |V^3|e^ 



;e+iAi 



IV^.I = ^j ri 1/^ (1 + ta^h (^^(7^,(x, y) - 0))) , 



A = -— (sin 0, — cos 0) , 



(Bl) 



where j = 1, 2, 3 and Uj is the ground state value of each 
superfluid density. The parameter ^j gives the core size while 
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6 and IZ are 



A^. 



^{x,y) = ^Gi{x,y), 



The initial position of a vortex is given by (x^, i/i). The func- 
tions Aab = ^b — ^a can be used to initiate a domain wall, 
when the ground state exhibits /7(1) x Z2 symmetry. 



Qi{x,y) = tan ^ 



X — Xj 



N^ 



n{x,y) = ^Uiix^y). 



^i{x,y) = y^{x - XiY + (y - yiY . 
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